Accepted for publication in the Astrophysical Journal 



Biases on initial mass function determinations. II. Real multiple 
systems and chance superpositions^ 

J. Mafz Apellaniz^'^ 

Instituto de Astrofisica de Andalucia-CSIC, Camino bajo de Huetor 50, 18008 Granada, 

Spain 

ABSTRACT 

When calculating stellar initial mass functions (IMFs) for young clusters, one 
has to take into account that most massive stars are born in multiple systems 
and that most IMFs are derived from data that cannot resolve such systems. It 
is also common to measure IMFs for clusters that are located at distances where 
multiple chance superpositions between members are expected to happen. In 
this article I model the consequences of both of those phenomena, real multiple 
systems and chance superpositions, on the observed color-magnitude diagrams 
and the IMFs derived from them. Using numerical experiments I quantify their 
influence on the IMF slope for massive stars and on the generation of systems 
with apparent masses above the stellar upper mass limit. The results in this 
paper can be used to correct for the biases induced by real and chance- alignment 
multiple systems when the effects are small and to identify when they are so large 
that most information about the IMF in the observed color-magnitude diagram 
is lost. Real multiple systems affect the observed or apparent massive-star MF 
slope little but can create a significant population of apparently ultramassive 
stars. Chance superpositions produce only small biases when the number of 
superimposed stars is low but, once a certain number threshold is reached, they 
can affect both the observed slope and the apparent stellar upper mass limit. 



^This article is partially based on observations made with the NASA/ESA Hubble Space Telescope (HST), 
some of them associated with GO program 10602 and the rest gathered from the archive, obtained at the 
Space Telescope Science Institute, which is operated by the Association of Universities for Research in 
Astronomy, Inc., under NASA contract NAS 5-26555. 

^e-mail contact: jmaiz@iaa.es. 

■^Ramon y Cajal fellow. 
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In the second part of the paper, I apply the experiments to two well known 
massive young clusters in the Local Group, NGC 3603 and R136. In both cases 
I show that the observed population of stars with masses above 120 Mq can 
be explained by the effects of unresolved objects, mostly real multiple systems 
for NGC 3603 and a combination of real and chance- alignment multiple systems 
for R136. Therefore, the case for the reality of a stellar upper mass limit at 
solar or near-solar metallicities is strengthened, with a possible value even lower 
than 150 Mq. An IMF slope somewhat flatter than Salpeter or Kroupa with 7 
between —1.6 and —2.0 is derived for the central region of NGC 3603, with a 
significant contribution to the uncertainty arising from the imprecise knowledge 
of the distance to the cluster. The IMF at the very center of R136 cannot be 
measured with the currently available data but the situation could change with 
new Hubble Space Telescope (HST) observations. 

Subject headings: binaries: general — binaries: visual — methods: numerical — 
methods: statistical — stars: luminosity function, mass function — open clusters 
and associations: individual (NGC 3603, 30 Doradus) 



1. Introduction 

This paper is the second one of a series where we explore the effec ts of different biases on 



the d etermination of the stellar initial mass function (IMF) . In paper I (IMaiz Apellaniz fc Ubeda 



20051 ) we analyzed the numerical biases induced by using bins of equal width when fitting 
power-laws to binned data (an effect that is more general than its application to the calcu- 
lation of mass functions). Those biases can be eliminated in several ways, of which a simple 
one is by grouping the data in equal-number bins (as opposed to equal- width bins). In this 
second paper I explore the effect of unresolved multiple systems, either physical or chance 
alignments, especially for the high-mass end of the IMF. In a future paper we will analyze 
the effect of random uncertainties in the mass determinations on the computed initial mass 
functioiil]. 



^Most of the analysis in these papers is also applicable to the calculation of the present-day mass function 
(PDMF). 



- 3 - 



1.1. Binaries and the IMF 



Why worry about binaries and higher-order multiple systems when calculating the IMF 
of a young stellar cluster such as NGC 3603 or R136? The quick answer is because one 
expects most of the observed point sources in clusters beyond a few kpc to be multiple. 
Therefore, deriving their masses by assuming they are single stars should introduce biases 
in the determination of the stellar mass functions of the clusters. 



Mason et al.l (119981 ) studied nearby O stars from a sample that included field objects 
and stars in low- and intermediate-mass clusters and associations. They analyzed the ob- 
served multiple fractions as a function of environment, including both spectroscopic and 
visual/speckle binaries and discovered that at least 72% of the O stars in clusters and asso- 
ciations had detected companions. Furthermore, the period distribution had a clear bimodal 
structure, with one peak around 10 days and another one around 10^ years. They argued 
that such an effect was quite likely to be instrumental because binaries with periods P near 
1-1000 yea rs are hard to detect and because the period distribution for B stars is quite fiat|^. 
Therefore, iMason et al.l (119981 ) predicted that when better spatial resolution were achieved 
and more intensive radial-velocity variation surveys were undertaken in the future, the gap 
between the two peaks might be filled and the true multiple fraction for massive stars in 
clusters might be found to be close to unity. Another interesting result of that paper was 
that the mass-ratio distribution for O binaries was fiat or nearly so, i.e. massive stars tend to 
have companions that are more massive on average than what would be expected from a ran- 
dom sampling of the IMF. The shape of the secondary mass distribution function is relevant 
because different slopes should af fect the observed or app arent mass function derived from 
unresolved data in different ways. ISagar fc Richtlerl (Il99ll ) calculated that the apparent MF 
slope for intermediate-mass stars in the random pairing case becomes significantly flatter 
when adding unresolved binaries but the effect is likely to be different in the fiat mass-ratio 



case. 



In the ten years a fter iMason et al.l (119981 ) other works have validated their conclusions. 
Garcia fc MermilliodI (120011 ) measured that 11 out of 14 0-systems in the cluster NGC 6231 
are spectroscopic binaries. This means that at least 79% of its 0/WR systems are mul- 
tiple and at least 88% of their 0/WR (and some early B) stars are in multiple systems. 
Their numbers are also similar fo r early -B stars. Also, new studies are indeed filling the 
Mas on period gap: iGaraen et al.l (120071 ) reported 11 new massive spectroscopic binaries 
and iMafz Apellaniz et al.l J2008b[ l cite two new visual systems in the solar neighborhood 



^Such a distribution, ex pressed in te rms of separation r instead of P, is close to /(logr) oc (logr)" and is 
known as Opik's law, after lOpikI (jl924f ). 
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see also iNelan et al.l I2004J '). With respect to the secondary mass distribution function, 



Kouwenhoven et aLl (|2005| ) studied intermediate-mass stars in the Sco 0B2 association and 
found that they could exclude random pairing in binaries, since the more massive stars 
cle arly tend to associate with other massive objects. Another recent study of Cyg 0B2 
by iKobulnicky fc Fryerl (120071 ) also finds a large (>80%) binary fraction for massive stars. 
Those authors favor a secondary mass distribution function that is either (a) nearly fiat 
or (b) a combination of a 60% Kroupa-like distribution and a 40% "twin" (M2 > 0.95Mi) 
component. 

But what about more massive clusters? Is it possible that the multiple fraction there 
is significant l y low er? Studying distant objects is more complicated but a recent survey by 
Evans et all J2006h found that NGC 346 (in the SMC), NGC 2004, and Nil (in the LMC) 
have a minimum spectroscopic binary fraction for O- and early B-type stars of 23 to 36%. 
Considering that their program did not allow them to detect binaries with periods of more 
than several tens of days and that at the distance of the Magellanic Clouds it is very hard 
to resolve visual binaries, the real binary fraction should be considerably higher. Therefore, 
it is safe to conclude that the effect of unresolved massive binaries in the observed IMF is 
important and needs to be analyzed. Such an analysis is the first goal of this paper. 



1.2. The stellar upper mass limit 



From the point of view of the origin of the IMF, it has become apparent in the last 
decade that massive stars must form in a different way than low-mass objects because of 
the short time scales involved and the limiting effect of radiat ion pressure on the mass 
accretion rate (IBally fc Zinneckerl l2005l : IZinnecker fc Yorkd 120071 ). There are currently two 
viable alternatives to form a massive star: by accretion through a disk, either in a relatively 
isolated environment or in a more crowded one where competition for mass among different 
protostars may be important, and by stellar collisions followed by mergers. The second 
alternative can only be relevant when the stellar densities are very high, such as at the core of 
a dense and massive cluster, but cannot be significant in relatively low-density environments 
such as OB associations. 



If stellar mergers are indeed impo rtant at the core of den s e clus ters, there should be large 
consequences on their stellar IMFs. iPortegies Zwart et al.l (119991 ) predicted that runaway 
collisions in such cores will lead to stars with masses above 100 Mq in a time scale of one 
or a few million years, which is comparable to the expected life times of those stars. The 
process may keep going on to the point of producing objects above 1000 M© that would 
end up as intermediate-mass black holes (IPortegies Zwart et al.ll2004j ). Therefore, massive 
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dense clusters should show an stellar upper mass limit mmax significa ntly higher than tho se of 
lower density or niass. Can such stellar behemoths exist? Recently, iBelkus et al.l (120071 ) and 
Yungelson et al.l (120081 ) have attempted modeling of solar- metallicity stars up to 1000 Mq, 
which would be formed by mergers. However, such modeling encounters problems with 
the determination of the Eddington limit, its associated instabilities and interaction with 
rotation, and with the large extrapolation required for mass-loss rates (which are already 
uncertain at the present time even for run-of-the-mill main-sequence O stars). Hence, theory 
has not given its final word on the possible existence of such ultramassive stars. 

What do observations tell us about those hypothetical ultramassive stars? Several pa- 



pers on R136 and other Local Group m assive young clusters (IWeidner fc Kroupall2004j : iFiger 



20051 : lOev fc Clarkell2005l : lKoenll2006f ) indicate that m^ax ~ 120-200 Mq when evolutionary 
tracks are used to obtain stellar masse^. Some luminous blue variables (LBVs) such as 
1] Car could have slightly higher evolutionary masses but LBVs are notoriously difficult to 
observe in detail due to the existence of circumstellar material. With respect to the more 
accurate masses measured from the orbital motion of the stars (keplerian masses), the high- 
„„uK„v,„j „„„„„„„„„j ^„ Asro on„ „ u;„„„y ^j^j^ 33 g ± 5.0 Mq and 82.0 ± 5.0 Mq 



(Bonanos et al. 


2004 


, see also 


Rauw et al. 


2004) 



that observations do not favor the existence of ultramassive stars at near-solar metallicity. 

Unresolved multiple systems can also play a role in the identification of ultramassive 
stars by making e.g. a binary made out of two stars close to the stellar upper mass limit 
appear as a single object w ith an evolutionary mass well above that limit. For example, 
Maiz Apellaniz et al.l (120071 ) found that Pismis 24-1, formerly an ultramassive star candidate, 
is in reality composed of three stars with masses below 120 Mq. It is also possible that the 
observed rrimax measured from evolutionary masses is affected not only by real (bound) 
multiple systems but also by chance alignments in the crowded environments where some of 
these objects are found. Analyzing the possibihty that some of the objects observed close 
to rrimax may be blended sources composed of stars of significantly lower mass is the second 
goal of this paper. 



2. Experiments 

In order to test the effect of unresolved binaries and blending on the massive-star IMF 
slope and the stellar upper mass limit I have designed three numerical experiments. The 
first one analyzes the effect of real binaries, the second one that of chance superpositions. 



■^Those values are called evolutionary masses and refer to the initial mass of the star. 
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and the third one combines the two effects. The purpose of the experiments is to analyze the 
effects on a simple test IMF, provide a way to correct them, and describe the tools needed 
to extend the technique to other circumstances. 

I design the simple case by assuming a well-sampled single stellar population with an 
age of 1 million years and solar metallicity. The choices for the age and metallicity are given 
by our desire to test our results with an analysis of NGC 3603 and R136 (note that the latter 
is currently thought to be of somewhat lower metallicity and may be slightly older) and by 
the simplicity of not having objects clearly evolved off the main sequence or even collapsed. 
Since we are only interested in the behavior of the top part of the IMF, for simplicity I assume 
that all the stars have already reached the main sequence, even though such an assumption 
should not valid for the low-mass stars. The Teg, logg, and log L information was obtained 



from the Geneva evolutionary tracks and isochrones database by_ Lejeune &: Schaererl (120011 ) 



e by 
tiorE 



selecting the case with stan dard mass-loss rates and no rotatioiu. For masses below 0.8 M© 
I use the My calibration of iKroupa et al.l (119931 ) to extend the isochrone until 0.1 M©. For 



the value of mmax, I select the highest value available in the Geneva tracks, 120 Mq, which 
is on the low end of the possible range (for mmax) for solar and near-solar metallicities. 
The 1-millio n-year isochrone was inserted into the new evolutionary synthesis module of 
CHORIZOS ( Maiz Apellaniz 2004) to calculate the Johnson U and V absolute magnitude^ 
as a function of mass usin g different synthetic spectra l energy distributions (SEDs) as a 



and Lejeune ( iLejeune et al 



function of T^ f. TLUSTY (Lanz fc Hubenvl l2003l. |2007|), Kurucz (jCastelh fc Kuruczl l2003l ) 



Il998h for the high, int ermediate, and l ow te mperature ranges 



respectively. The Joh nson zero points derived by iMaiz ApellanizI (120061 ) and modified by 
Maiz Apellanig (120071 ) were used. The final results are two tabulated, finely-gridded functions 
that give V and U — V as a. function of mass m between 0.1 and 120 M©. The two functions 
were then extrapolated in logm up to m = 1000 Mq (logm = 3). Note that V{m) is a strict 
monotonically decreasing function over the full range of m = 0.1 — 1000 Mq (this is true for 
the chosen isochrone but it is not true in general). Therefore, it is possible to invert it to 
produce m{V), which gives the true mas^ that corresponds to a star of absolute magnitude 
V. 



"^In the last years much work has been done on new evolutionary tracks and isochrones with rotation and 
possibly-reduced mass-loss rates due to clumping corrections. The present work can be easily adapted to 
those other conditions but the effects on our results are expected to be relatively minor. 

^For the sake of notation simplicity and since at this point I am not discussing distance effects, I write 
Mjj and My simply as U and V, respectively. 

^Throughout this paper m is considered to be the initial stellar mass, not the mass when the star is 1 
million years old. 
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Once V{m) and U{m) — V{m) as a function of mass have been generated, I calculate 
the corresponding bidimensional single-star color-magnitude density function gis{U — V, V) 
by applying one of two cases: [a] a Kroupa IMF, which can be described as a continuous 
function comprised of two power law segments dn/dm = f{m) oc m'^, with 7 = —1.3 for 
m = 0.1 - 0.5 Mq and 7 = -2.3 for m = 0.5 - 120 Mq; [b] a flatter, top-heavy IMF, similar 
to Kroupa but with 7 = -1.3 for m = 0.1 - 0.5 M© and 7 = -2.0 for m = 0.5 - 120 Mq. 
9isiU — V, V) can be collapsed into a one- dimensional magnitude functioiill G'is(V^) without 
losing any information because for our isochrone V{m) is a strictly monotonic function. Since 
both gis{U — V,V) and GisiV) are functions defined in a fine (2-D and 1-D, respectively) 
grid, special care has to be taken to avoid introducing numerical biases and noise due to 
sampling and interpolation. 

Finally, in order to reproduce the blending problems associated with the strategy com- 
monly used to derive evolutionary masses from luminosity functions, I combine individual 
stars into blended objects according to different rules. This is done by adding their U- 
and y-band luminosities and computing the resulting V and U — V values. From there, 
I generate: [a] a new bidimensional color-magnitude density function gx{U — V,V) from 
gis{U — V,V), where x specifies the combination rule, and [b] the apparent mass function 
fa,x{^a) calculated by collapsing gx{U — V, V) into G^iV) and then converting the observed 
absolute magnitudes into apparent masses. The last step is done by applying m{V), which is 
a relationship derived for individual stars that is not correct for blended objects. Therefore, 
the derived apparent masses are not the real ones. In other words, I am deriving an 
apparent mass function (AMF) by [a] assuming that all the observed objects are single stars 
(even though they are not) and [b] using only the absolute magnitudes (but not their colors) 
to obtain the masses. The first assumption is applied in most IMF studies. The second 
assumption should be valid as long as all the stars are at the same distance, extinction is 
well known, contamination by field objects or stars of different ages is negligible, and m{V) 
is a strictly monotonic function. Note, however, that for ages older than 2.5 million years 
(when mlV) is no longer monotonic) the same type of analysis can still be done measuring 
masses by minimizing the distance to a given isochrone. 

Note also that I am not using spectral types for the calculation of the IMF. For massive 
stars, the main advantages of using spectral types to calculate IMFs are a better determi- 
nation of the cluster age, a direct knowledge of the intrinsic color of the star (hence, usually 
a better extinction measurement), and an ability to derive additional luminosity/distance 
information. Those are all possible sources of biases, so spectroscopy indeed represents an 



^In this paper I use magnitude functions instead of luminosity functions; the conversion between one and 
the other is straightforward. 
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improvement over a pure photometric analysis. However, as previously stated, for the exper- 
iments in this paper I assume that age, extinction, and distance are all well constrained by 
external information, so those advantages of spectroscopy are not relevant for our discussion. 
The additional advantage that spectroscopy has over a simple color-magnitude analysis when 
blended sources are present is the possibility of identifying them by their composite spectra 
(something that can also be done in some cases with multi- filter photometry). Note, how- 
ever, that many spectra of blended sources cannot be distinguished from those of single ones 
for three possible reasons: [a] need of very high S/N data to detect the contribution from 
a dim companion, [b] absence of time-resolved data to detect spectroscopic binaries, and [c] 
similarity between the two blended spectral^. In summary, spectroscopy can be and in most 
cases is a helpful aid to the elimination of biases in IMF determinations but it also has its 
costs and limitations in detecting blended sources. This paper should not be interpreted as 
neglecting the importance of spectroscopy for the study of the massive-star IMF but rather 
as a study of what can be done when spectroscopy is not available. 



2.1. Experiment 1: Real binaries 

For the first experiment I attempt to reproduce the effect of unresolved multiple systems 
which, for simplicity, I assume to be always binaries. The combination rules applied in this 
case are: 



1. Each star is blended with one and only one star. 

2. For a blended system with mi > m2, m2 is randomly selected from a flat distribution 
in mass between 0.1 M© (the lower mass limit) and mi. The pairing is done while 
keeping the real, single- mass stellar IMF as Kroupa (case a) or top-heavy (case b). 



This experiment reproduces the physical situation in which all the stars in a cluster are 
binaries with a fiat mass-ratio distribution. As previously discussed, this is a rather accurate 
approximation for systems with at lea st one massive s t ar bu t is likely to be more imprecise 
for lower-mass objects. For example, iReipurth et al.l (120071 ) o bserve only an 8.8% binary 



fraction among low-mass stars in the Orion Nebula Cluster and ISoUima et al.l (120071 ) detect 
a binary fraction of ~ 40% for globular clusters with ages less that 10^° years (of course, the 



^For example, the spectra of the two components of WR 20a (|Bonanos et al.l 120041 ) are almost indistin- 
guishable and if it were not for the large velocity variations induced by their proximity it would be almost 
impossible to identify the spectrum as a composite. 
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Flat mass-ratio binaries Kroupa 




Fig. 1. — Color- magnitude density function gihiU — V, V) for the Kroupa case of experiment 
1 (real binaries) shown as a Hess diagram. The black line shows the position of the 1 million 
year isochrone. The function scaling is logarithmic and the normalization is arbitrary. See 
the electronic version of the journal for a color version of this figure. 
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situation in globular clusters is complicated by dynamical evolution, as the lower fraction 
detected by the same authors for older clusters show). However, since our goal is to study the 
top part of the resulting AMF, the validity of the approximation for low-mass objects should 
not concern us much. The other assumption in this experiment is that all massive binaries 
are unresolved. Under what circumstances is that reasonable? As previously described, in 
the solar neighborhood (< 2 kpc) ~ 1/3 of the massive systems in clusters/associations 
are spectroscopic binaries, another ~ 1/3 are visual/interferometric multiple systems, and 
the remaining ~ 1/3 could be undetected multiples (with periods approximately between 
those of the detected spectroscopic and visual/interferometric systems )0. This means that 
even in the solar neighborhood ~ 2/3 of the massive binary systems would appear as point 
sources in a ground-based imaging study with standard imaging techniques. Using better 
spatial resolution (e.g. HST or non-standard techniques such as adaptive optics or lucky 
imaging) buys you a larger fraction of resolved systems in the solar neighborhood but once 
you reach several kpc, even that is not enough. For example, ACS/HRC (the imager with 
the best spatial resolution on HST) can still resolve objects with a magnitude difference 
Am f» 1 mag and a separation of 50 mas but that is close to the limit of its capabilities 



Mafz Apellaniz et al.ll2008bl ). At a distance of 7 kpc, 50 mas corres ponds to 350 



(see e.g. 

AU, w hich is the typical value for the separation in the visual 0-systems of iMason et al. 



( 119981 ) ■ implying that at that distance ACS/HRC would resolve only 20-30% of the massive 
binary systems (if their properties are similar to those in the solar neighborhood) . The LMC 
is hopeless, even with HST: 50 mas corresponds to 2500 AU, so only very-large separation 
systems are resolved. Therefore, I can conclude that this experiment represents a reasonable 
approximation to the conditions in many IMF studies of young clusters and associations 
beyond the solar neighborhood. 

The resulting color-magnitude density function, gih{U — V,V) is shown in Fig. [1] for 
the Kroupa IMF case (the top-heavy equivalent function is not shown but is visually very 
similar). gih{U — V,V) is zero nearly everywhere in the color- magnitude diagram except 
in the narrow band between the original isochrone and a near-parallel line that is located 
approximately 0.75 magnitudes above it (the location of the equal-mass binaries). The band 
between the two lines is filled with the non-equal-mass binaries. At the low-mass extreme 
the band becomes narrower and detached from the isochrone because the addition of any 
companion changes the magnitude (and possibly the color) considerably. That effect does 
not take place for higher masses because the addition of a 0.1 Mq star is undetectable in 
color-magnitude with the scale of Fig. [1] At the high-mass end an obvious extension beyond 



9 



It should be pointed out that th e ov erlap between spect roscopic and visual/interferometric systems is 



very small, see lBovaiian et al.l (|2007t ) and lNorth et al.l (|2007r) for examples 
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Fig. 2. — (Real) IMF and AMF for experiment 1 (real binaries). The top panel shows the 
Kroupa case and the bottom panel the top-heavy case. The mass is expressed in solar units. 
The IMF is normahzed to 1.0 and the AMF to 0.5. 
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^max is observed. The extension is nearly vertical because massive stars of this age have 
near-degenerate U — V colors. The apparent upper mass limit, mmax,ib is well defined and 
has a value of 182 Mq. 

The calculated AMFs for the Kroupa and top-heavy cases are shown in Fig. [2] along 
with the original IMFs. The behavior in both cases is similar. For low masses, the apparent 
lower limit is slightly shifted towards higher values (corresponding to a 0.1 Mq + 0.1 Mq 
binary, since 0.1 Mq is the lower mass limit of our experiments) and several peaks in the 
AMF can be seen below 0.7 Mq due to the fine structure of m(y)0. For intermediate and 
high masses, the AMF runs nearly parallel to the IMF at a value roughly 1/2 of the IMF (the 
AMF includes one half of the systems of the IMF). In order to quantify the change in slope 
between the IMF and the binary AMF, A7, I have fitted a power law to a number of mass 
ranges using the technique in Paper I. Results are plotted in Fig. [31 For miower between 1 and 
16 Mq and mupper 7^ 120 Mq there is a steepening of the MF with A7 between —0.2 and 0.0 
while for larger masses in most cases there a slight fiattening with A7 between 0.0 and 0.2. 
A7 depends only weakly on mupper and the IMF slope. Note that some of the small-scale 
structure observed in Fig. [3] at the 0.02 level may be of numerical origin. I conclude that for 
most cases the existence of unresolved binaries has only a small effect on the massive-star 
IMF slope. 

The largest difference between the IMF and the AMF happens beyond rumax, where a 
tail is formed in the AMF: these are the apparently ultramassive stars (AUMSs). In Tabled] 
I quantify the effect by giving the apparent number ratios between the objects above mmax 
and the objects in different ranges below it. There is a strong dependence on the IMF slope, 
since a top-heavy IMF allows a significantly larger number of objects to become AUMSs. 

Figure [T] reveals another interesting issue which is relevant to IMF studies of Galac- 
tic clusters, where distance is one of the unknowns that has to be solved for. O stars are 
notoriously difficult to use as yardsticks because their near-degenerate optical colors place 
them along near-vertical liii es in a color-magnitude diagram (see e.g. Fig. [1] but see also 



Maiz Apellaniz fc Sotall2008l for recent progress on improving our knowledge of 0-stars col- 
ors). For that reason, main-sequence B stars (which have intrinsic U — V in the approximate 
range from —1.3 to —0.1) are more commonly used to derive distances with main-sequence 
fitting techniques. However, if the assumption which underlies our experiment is relevant 
(that most O and B stars are unresolved binaries with fiat mass-ratio distributions), then it is 



^"The conversion from a magnitude or luminosity function to a IMF involves dm/dV, whose numerical 
behavior may be odd when using coarsely-gridded tabular information to derive it. The observed structures 
may or may not be real but since our concern is with the high-mass end of the AMF, we will not discuss 
them any further. 
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lower ^ Sol' 



Fig. 3. — Slope change introduced by flat mass-ratio binaries (first experiment) as a function 
of the lower limit for the mass range (miower) for the two IMF cases and for different values 

). In each case the values between 1 Mq and 



m 



upper 



of the upper limit for the mass range 
0.60mupper are plotted. See the electronic version of the journal for a color version of this 
figure. 
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easy to see that a direct application of main-sequence fitting to aU — V vs. V diagram would 
lead to an overestimation of the distance, because for a fixed value of U — V there should be 
more stars close to the equal-mass binary sequence than to the main sequence itself. Note 
that this effect should be significantly smaller for older clusters be cause the binary fra ction is 
lower there, as it is indeed observed for globular clusters (see e.g. ISollima et al.ll2007l ). Spec- 
troscopy can partially resolve the problem by detecting blended spectra or radial velocity 
variations but from Fig. [1] we can see that it should be also possible to correct for the effect 
by considering that there is an intrinsic spread in the observed color-magnitude diagram 
due to binarity (besides the ones due to differential extinction and foreground/background 
contamination). 



2.2. Experiment 2: Chance superpositions of single stars 

In the second experiment I model the effect of chance superpositions between single 
(not binary) stars within the observed clusters. Such a simulation introduces additional 
complications that are detector- and observation-dependent because now we have to deal 
with a continuous spatial distribution of stars along a 2-D detector. In other words, the 
observed luminosity function depends not only on the intrinsic properties of the cluster (e.g. 
its internal structure) but also on the the relative pixel and point-spread function (PSF) sizes, 
Poisson statistics, and the ability of the detector and the finding algorithm to differentiate 
between single- and multiple- point sources of diverse Am. Given those considerations, the 
goal in this case is to obtain a toy model that allows us to approximately quantify the effect 
for a given region of a cluster in order to test the magnitude of the effect. 

The combination rules in this case are: 



1. The light from a given star falls in a single pixel (the fat pixel approximation) . 

2. The initial cluster is built by generating a fine grid inU—V and V and assigning to each 
point a value of gis{U — V, V) integrated over the size of the bin. Since gis{U — V, V) 
is zero outside the 1-million year isochrone, most points in the grid have no stars. The 
initial cluster is then defined by the non-zero values of gis{U — V,V). As with the 
previous experiment, we consider the two cases of a Kroupa and a top-heavy IMF. 

3- g2s{U — V,V) is a color-magnitude function that describes a cluster with two single 
stars in each pixel which are randomly paired. I generate it by convolving gis{U — V, V) 
with itself and assigning the result of each pair to the closest value in the {U — V, V) 
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CTl Note that g2s{U — V, V) does not take into consideration Poisson statistics for 
the stellar distribution among pixels but instead considers that each pixel contains 
exactly two stars. 

4. g4s{U — V,V) is generated by convolving g2s{U — V,V) with itself, and in general, 
g2"s{U — V, V) is generated by convolving g2^-^siU — V, V) with itself. In this paper I 
generate the color-magnitude functions up to n = 11 i.e. 5'2048s(f^ — V,V)- 

Those rules should be able to simulate the general behavior of the observed color- 
magnitude diagram as the intrinsic stellar density or the cluster distance increases. However, 
there are three caveats to its direct application. The first one is that a comparison with a real 
case should always consider the possible additional effect of Poisson fluctuations, especially 
when the number of superimposed stars A^^sup = 2" is small. One way to do this if the expected 
average number of stars per pixel, iVpx, falls between 2"~^ and 2" would be to consider both 
of those possible. The second caveat is that the fat pixel approximation is not valid 

for a well-sampled detector (PSF size several times larger than the pixel size), so an effective 
pixel size has to be defined by calculating the area of influence of a star i.e. the region 
where the presence of a star does not allow a second star to be detected but instead the 
fluxes from both are merged. Such an effective pixel size is in reality Am-dependent so an 
average effective pixel size should be calculated. Also, the precise flux derived depends on 
the photometric technique (PSF fitting or aperture photometry) and the algorithm details. 
Once again, I stress that the purpose of this experiment is not to find the exact correction 
required to account for chance superpositions but rather to develop a toy model that would 
allow the observer to roughly quantify the effect. The third caveat is that all stars are 
assumed to be single. This last caveat will be addressed by the third experiment. 

Six of the resulting full color-magnitude density functions, g2^siU — V,V), are shown 
in Fig. m for n = 1,4,7 and the two cases (Kroupa and top heavy). The top left corner 
(where massive stars are located) of all the functions is shown in Figs. [5] (Kroupa) and [6] 
(top heavy). As was the case in experiment 1, g2^^s{U — V, V) is zero nearly everywhere in 
the color-magnitude diagram except in a narrow band. The band widens up until n = 8 and 
then narrows down from that point onwards. At the same time, the lower end of the band 
starts climbing up the vertical coordinate starting at n = 1. The top end starts climbing up 
at a lower pace but then picks up speed around n = 8. The progression of both ends is faster 
for the top-heavy case than for the Kroupa case. For the final n = 11 value the resulting 
Hess-diagram has been compressed to a well-defined narrow peak in the U — V vs. V color 



As previously noted, special care was taken to avoid roundoff errors. 
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Fig. 4. — Color-magnitude density functions g2s{U — V,V) (left), gi^siU — V,V) (center), 
and gi2Ss{U — V, V) (right) for the Kroupa (top row) and top-heavy (bottom row) cases of 
experiment 2 (chance superpositions of single stars) shown as Hess diagrams. The thick black 
line shows the position of the 1 million year isochrone. The function scaling is logarithmic 



and the normalization is arbitrary. See the author's web site http://www.iaa.es/~jmaiz 

for animated gifs showing all the color-magnitude functions g2"s{U — V,V) from n 
n = 11. See the electronic version of the journal for a color version of this figure. 
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Fig. 5. — Top left comer of the color-magnitude density functions g2-^s{U — V,V) with 
n = 0, 11 for the Kroupa case of experiment 2 (chance superpositions of single stars) shown 
as Hess diagrams. The n = case has been smoothed in order to allow it to show some 
extent perpendicular to the 1 million year isochrone. The thick black line shows the position 
of the 1 million year isochrone and the thin lines the evolutionary tracks between and 2 
million years for the initial masses of 25, 40, 60, 85, and 120 Mq. The function scaling is 
logarithmic and the normalization is arbitrary. See the electronic version of the journal for 
a color version of this figure. 
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Fig. 6. — Top left comer of the color-magnitude density functions g2-^s{U — V,V) with 
n = 0, 11 for the top-heavy case of experiment 2 (chance superpositions of single stars) 
shown as Hess diagrams. The n = case has been smoothed in order to allow it to show 
some extent perpendicular to the 1 million year isochrone. The thick black line shows the 
position of the 1 million year isochrone and the thin lines the evolutionary tracks between 
and 2 million years for the initial masses of 25, 40, 60, 85, and 120 Mq. The function scaling 
is logarithmic and the normalization is arbitrary. See the electronic version of the journal 
for a color version of this figure. 
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diagram: the distribution is completely contained in the last plot of the top-heavy case of 
Fig. [6] and almost so for the Kroupa case of Fig. [51 

FiguresHlEl andOindicate that, as crowding becomes more of an issue, we progressively 
lose the ability to detect low-mass stars. Once a certain threshold has been attained, high- 
mass stars also become significantly affected by crowding. However, one should be careful 
in overinterpreting those Hess diagrams as color-magnitude probability distributions for the 
observed population of points source in a a given cluster because of the fat pixel assumption. 
The lower end of the observed band represents in reality those pixels which detection algo- 
rithms would simply classify as background. Therefore, a more realistic interpretation of the 
Hess diagrams would be that the observed point source population measured by assuming 
a zero background (or calculating it outside the cluster boundaries) is extracted from the 
top part of the 2-D function given by g2^s{U — V,V) and that the lower part represents 
the unobserved part of the population. The limit between the two regions is defined by the 
characteristics of the detector and the observations. 

An interesting alternative interpretation of Figs. HI and E] is that they represent the 
probability distribution for a cluster with 2" stars, with the extension in U — V and V 
indicating the result of stochastic effects in drawing the given number of stars from an IMF. 
I will not explore th i s inte rpretation any further in this article; the reader is referred to 



Cervino fc Luridianal (120061 ) for more details. 



It is also useful to compare the results from gi^ (first experiment) and g2s (second 
experiments)!^. Both functions represent the combination of two stars but with an important 
difference: in the first case the pairing assumes a fiat mass-ratio distribution while in the 
second case the pairing is random. The different combination rules manifest themselves in 
the slope that is obtained when plotting ^fib or g2s as a function of V for a fixed U — V. 
For most values of U — V in the first case the function increases towards lower values of V 
(higher luminosities) while in the second one the behavior is the opposite. The reason is that 
for random pairing, once an mi is drawn in the middle or upper main sequence, chances are 
that m2 will be significantly lower. On the other hand, for a flat mass-ratio distribution one 
gets a larger contribution from stars with 1712 only slightly lower than mi, thus tilting the 
balance at a flxed U — V towards systems with near-equal mass ratios. 

The calculated AMFs for the Kroupa and top-heavy cases are shown in Fig. [7] along with 
the original IMFs. The overall behavior in both cases again being similar: as the number 
of stars increases, the left extreme of the AMF moves towards the right in approximately 
equal amounts in logm for each step in n. The right extreme has a different behavior: the 



^^Note that the function scahng in Figs. [T] and U is different. 



-20- 



Table 1. Apparent ratios between the number of objects above 120 Mq and the number in 
different mass ranges for the two cases of the first experiment. 



Mass range 


Kroupa 


Top-heavy 


Me 


real IMF 


real IMF 


All 


7.873 • 10--'' 


3.639 • 10-^ 


1 - 120 


6.661 • 10-^ 


2.246 ■ 10-3 


3 - 120 


2.835 • 10-3 


7.142 • 10-3 


5-120 


5.941 ■ 10-=^ 


1.294 • 10-2 


8-120 


1.164 ■ 10-2 


2.225 • 10-2 


15 - 120 


2.921 ■ 10-2 


4.726 • 10-2 


25 - 120 


6.206 ■ 10-2 


8.895 • 10-2 


40 - 120 


1.284 • 10-1 


1.662 ■ 10-1 


60 - 120 


2.656 ■ 10-1 


3.163 • 10-1 



Table 2. Apparent ratios between the number of objects above 120 Mq and the total 
number of objects for the second (single stars) and third (binaries) experiments. 



#of 


Kroupa 


Top-heavy 


Kroupa 




Top-heavy 


objects 


sing 


les 


sing 


les 


binaries 




binaries 


1 


0.000 • 


10+0 


0.000 - 


10+0 


7.873 - 10" 


-5 


3.639 - 10-^ 


2 


2.080 ■ 


10-^ 


3.331 • 


10-6 


1.913-10- 


-4 


8.708 • 10-^ 


4 


1.027- 


10-6 


1.582 - 


10-5 


4.458 - 10- 


-4 


1.988 • 10-3 


8 


4.829 ■ 


10-6 


7.101 - 


10-5 


1.033 - 10- 


-3 


4.552 • 10-3 


16 


2.485 ■ 


10-5 


3.552 - 


10-^ 


2.469 ■ 10- 


-3 


1.082 • 10-2 


32 


1.422 • 


10-^ 


1.977- 


10-3 


6.209 - 10" 


-3 


2.703 - 10-2 


64 


9.067- 


10-^ 


1.155 - 


10-2 


1.646 - 10- 


-2 


7.090 - 10-2 


128 


6.040 • 


10-3 


6.253 - 


10-2 


4.626 - 10- 


-2 


1.948 - 10-1 


256 


3.697- 


10-2 


2.790 - 


10-1 


1.385 - 10" 


-1 


5.045 - 10-1 


512 


1.907- 


10-1 


7.765 - 


10-1 


4.079 - 10" 


-1 


9.118-10-1 


1024 


6.667 ■ 


10-1 


9.979 ■ 


10-1 


8.684 - 10- 


-1 


9.998 • 10-1 


2048 


9.944 - 


10-1 


1.000 - 


10+0 


9.997 - 10- 


-1 


1.000 - 10+° 
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Fig. 7. — IMF and AMFs for experiment 2 (chance superpositions of single stars). The 
top panel shows the Kroupa case and the bottom panel the top-heavy case. The mass is 
expressed in solar units. The IMF is normalized to 1.0 and the AMFs to 2~". 
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Fig. 8. — Slope change introduced by chance superpositions of single stars (second experi- 
ment) as a function of the number of superimposed single stars for the two IMF cases and 
for three different mass ranges. See the electronic version of the journal for a color version 
of this figure. 
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population of AUMSs initially appears only as a weak tail and it is only when = 7 — 8 
that a significant number is built up. From then on, the right extreme of the AMF picks up 
speed and starts moving towards the right, though at a lower pace than the left extreme. 
As a consequence, the distribution narrows in logm. The narrowing process also affects the 
overall shape of the distribution, which for n = 11 is already quite similar to a log-normal 
(especially for the top-heavy case). 

Table [2] lists the fraction of objects with apparent masses above m^^^ = 120 Mq for 
both the Kroupa (second column) and top-heavy cases (third column) as a function of A'sup- 
The fraction is highly dependent on the IMF slope, with the top-heavy case showing values 
an order of magnitude above the Kroupa case until the point where the value becomes close 
to unity for the top-heavy case {n ^ 9). Such a dependence is not an unexpected behavior, 
because even a small change in 7 is sufficient to produce a large variation in the ratio of 
stars with masses around 100 M© to low-mass stars in the real IMF (compare e.g. the two 
panels in Fig. [2]). As also expected, flat-mass ratio binaries are significantly more efficient 
than chance superpositions in creating a population of AUMSs: ~ 32 superimposed stars 
are required to create the same number of those stars as those created by a single binary for 
the Kroupa case (the number is slightly lower for the top-heavy case). 

Figure [8] plots A7 the change in MF slope as a function of A'sup for three mass ranges. 
Three regimes can be distinguished. For low A^sup, the effect in A7 is small and in most 
cases implies a slight steepening of order 0.1. For intermediate values, A7 starts to increase 
(the mass function becomes flatter) but the effect remains small and should be correctable. 
Finally, for large Agup, A7 becomes so large that no realistic correction is possible. 



2.3. Experiment 3: Chance superpositions of binaries 

The third experiment is a combination of the previous two. I start with the color- 
magnitude generated in the flrst experiment for flat mass-ratio binaries, gih{U — V,V), and I 
convolve it with itself to derive g2h{U — V,V). g^hiU — V, V) is then generated by convolving 
g2s{U — V, V) with itself and the process is repeated until n = 11 i.e. 5'2048b(t^ This 
experiment should provide the most realistic representation of the observed color-magnitude 
density function and AMF of a distant cluster, always keeping in mind the flrst two caveats 
previously discussed. 

As I did for the second experiment, I show six of the resulting full color-magnitude 
density functions, g2^h{U — V, V) in Fig. [9] for n = 1,4,7 and the two cases (Kroupa and top 
heavy). The top left corner (where massive stars are located) of all the functions is shown in 
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Fig. 9. — Color-magnitude density functions g2hiU — V,V) (left), giQhiU — V,V) (center), 
and gi2Sh{U — V, V) (right) for the Kroupa (top row) and top-heavy (bottom row) cases of 
experiment 3 (chance superpositions of binaries) shown as Hess diagrams. The thick black 
line shows the position of the 1 million year isochrone. The function scaling is logarithmic 



and the normalization is arbitrary. See the author's web site http://www.iaa.es/~jmaiz 

for animated gifs showing all the color-magnitude functions g2^h{U — V, V) from n 
n = 11. See the electronic version of the journal for a color version of this figure. 
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Fig. 10. — Top left corner of the color-magnitude density functions g2"h{U — V, V) with 
n = 0, 11 for the Kroupa case of experiment 3 (chance superpositions of binaries) shown 
as Hess diagrams. The n = case has been smoothed in order to allow it to show some 
extent perpendicular to the 1 million year isochrone. The thick black line shows the position 
of the 1 million year isochrone and the thin lines the evolutionary tracks between and 2 
million years for the initial masses of 25, 40, 60, 85, and 120 Mq. The function scaling is 
logarithmic and the normalization is arbitrary. See the electronic version of the journal for 
a color version of this figure. 
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Fig. 11. — Top left corner of the color-magnitude density functions g2"h{U — V,V) with 
n = 0, 11 for the top-heavy case of experiment 3 (chance superpositions of binaries) shown 
as Hess diagrams. The n = case has been smoothed in order to allow it to show some 
extent perpendicular to the 1 million year isochrone. The thick black line shows the position 
of the 1 million year isochrone and the thin lines the evolutionary tracks between and 2 
million years for the initial masses of 25, 40, 60, 85, and 120 Mq. The function scaling is 
logarithmic and the normalization is arbitrary. See the electronic version of the journal for 
a color version of this figure. 
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Figs. [To] (Kroupa) and [TT] (top heavy). Those three figures are quite similar to the equivalent 
ones for the second experiment. The most prominent differences are: [a] the non-zero regions 
of g2"h{U — V,V) are more extended in their short direction (in the diagrams, from the 
lower left to the upper right) compared to the non-zero regions of g2"s{U — V, V) due to the 
additional widening introduced by binaries; and [b] the extension of g2"h{U — V, V) in its long 
direction (from the lower right to the upper left) is similar to that of g2"s{U — V, V) but it is 
slightly displaced towards lower magnitudes (higher luminosities) due to the additional factor 
of two in the number of superimposed objects. The latter effect is especially noticeable for 
large values of n (compare the lower right panels of Figs. [5] and [10] or those of Figs. [6] and [TT]) . 
Note, however, that a comparison of the two experiments assuming the same number of 
superimposed individual stars implies using g2^s{U — V,V) and 5'2"-ib(f^ ~ V^V)- If the 
comparison is done in that way one finds that for large n the top parts of the non-zero regions 
of the two functions end near a similar value (with the g2"-^h{U — V,V) one only slightly 
higher in luminosity) and that the low-luminosity tail of g2"~^h{U — V,V) is significantly 
more extended than that of g2"s{U — V^V). 

The calculated AMFs for the third experiment are shown in Fig. [121 with the overall 
behavior in both cases being similar to the equivalent ones for the second experiment. A 
more detailed comparison shows that for low values of n the behavior beyond mmax is quite 
different: both experiments start by developing a tail from their n = case but since each 
initial distribution has a very different apparent upper mass limit (for experiment 2 it is the 
real mmax while for experiment 3 it is mmax,ib) the fraction of stars above the limit starts 
from a much higher value for experiment 3 and grows at a significantly slower pace (see 
Table [2]). For large values of n the left side of g2^-^h{U — V, V) is quite similar to that of 
5'2"-is(f^ — V, V) while its right side is similar to g2^s{U — V,V). 

Figure [13] plots the change in MF slope A7 as a function of N^^-p for three mass ranges. 
The same three regimes as in Fig. [H] can be distinguished. For low A^supj the effect in A7 
is small and can be approximated by the value derived from the first experiment. For 
intermediate values, A7 increases moderately (the mass function becomes fiatter) and the 
effect remains correctable. Finally, for large A^^sup, A7 becomes very large and no realistic 
correction is possible. 

3. Sample applications 

In this section I apply the results of the experiments above to the IMFs derived from 
HST observations for two of the best studied massive young clusters in the Local Group, 
NGC 3603 and R136. In order to do that, two preliminary steps are required: 
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Fig. 12. — AMFs for experiment 3 (chance superpositions of binaries). The top panel shows 
the Kroupa case and the bottom panel the top-heavy case. The mass is expressed in solar 
units. All AMFs are normalized to 2~^""^^^. 
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Fig. 13. — Slope change introduced by chance superpositions of binaries (third experiment) 
as a function of the number of superimposed binaries for the two IMF cases and for three 
different mass ranges. See the electronic version of the journal for a color version of this 
figure. 
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1. Derive the radius of an effective pixel. As previously mentioned, the value depends 
not only on the detector properties but also o n Am. Indeed, ACS/HRC on HST can 
separate low Am objects only 2 pixels apart (IMaiz Apellaniz et al.l l2008al ) (effective 
radius of 1 pixel) but for large Am the effective radius can grow to 4-5 pixels (similar 
values apply to the PC of WFPC2). As a compromise, we adopt a value of 2.5 pixels. 

2. Determine the average number of stars per effective pixel, A^^epx, given by: 



epx 



d 



10 pc 



j^g-0.4(mv,cpx-A^V,l) 



In Eqn. [H c? is the distance to the object, myepx is the averaged V^-band extinction- 
corrected measured magnitude in a single effective pixel, and Myi is the absolute 
magnitude of a "1-star cluster" derived from evolutionary synthesis (calculated by 
deriving My for a large cluster with e.g. = 10^ stars and then adding 2.5 log A^). If 
needed, V^-based values can be substituted by those from other passband. 



3.1. NGC 3603 



NGC 3603 is the massive young cluster with the lowest extinction in the Galaxy. 
Moffat et al.l (119941 ) described it as a Galactic clone of R136 without the halo due to their 
similarities in age, stellar content, and central density. Note, however, that most of the 
mass of 30 Doradus is located outside R13 6, making it si g nifica ntly more massive than NGC 
3603. A s tudy with HST spec troscopy by iDrissen et al.l (Il995l ). recently extended from the 
ground by lMelena et al.l (120071 ) . indicates that NGC 3603 contains the highest concentration 
of well-classified 02/3 -|- WNha stars in the Galaxy. 



I show in Table 13.11 the literature data for the distance and extinction to the core 
of NGC 3603. The values for the extinction are quite homogeneous and here we adopt 
Ay = 4.50. Also, at the location of the central cluster there is little poi nt-to-point variation 
in the extinction ( Moffat 1983 : Melnick et al. 19891 : Melena et al. 2007 ). so the assumption 
of a constant foreground screen appears reasonabk If]. The disagreement on the distance 
to NGC 3603 is la rger. Different met hods and observations yield values between ^ 6 kpc 



and ~ 8 kpc (see iMelena et al.l 120071 for a discussion). In order to consider the different 



possibilities, I adopt three possible distances of 6, 7, and 8 kpc, which I refer to as the short, 
intermediate, and long distances, respectively. Regarding the age of the core of NGC 3603, 



Strong extinction variations exist at larger radii but those do not concern us here. 
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it was previously thought to be ~ 3 milhon years based on the presence of WR stars but 
now it is beheved that th ose objects are s till in their core hydrogen-burning stage (i.e. they 
are WNha stars, see e.g. ICrowtherl 120071) and the currently favored age is 1 million years 
Jstolte et ahlbooi ISung fc Besselllbooi ). which I adopt. 



NGC 3603 was observed with HST/ACS on 29 December 2005 under GO program 10602 
using both the WFC and HRC detectors and several filters. Here I use the HRC F550M data, 
combined with the literature information described above for the distance, extinction, and 
age, to study the blending effect created by real multiple systems and chance alignments on 
the NGC 3603 IMF. In a future paper we will use all the data to do a more thorough analysis 
that will include independent measurements of the distance and extinction. Four F550M 
exposures of 2 s each were obtained with a dithering pattern that covers the gap produced 
by the HRC occulting finger. The exposure time was chosen in order to avoid saturation 
of the brightest stars. MULTIDRIZZLE was used to generate a common geometrically- 
undistorted cosmic-ray-corrected drizzled (DRZ) frame, see Fig. [TH as well as to clean the 
individual geometrically-distorted flat-fielded (FLT) frames from cosmic rays and hot pixels. 

JMAPHOT, a crowded-field photometry package specifically designed for HST images 
was used for the data reduction. First, a search is performed in the DRZ frame to find all the 
sources above a certain S/N threshold. JMAPHOT finds the sources in order from bright 
to dim and, after each one is found, a model (the drizzled PSF for point sources, simple 
kernels for extended sources and possible leftover cosmic rays and hot pixels) is subtracted 
from the data and a correction applied to the weight map before the next source is searched 
for. This technique allows for the exclusion of false positives caused by complex PSFs with 
secondary peaks and for the detection of some of the otherwise false negatives close to bright 
sources. Once a source list is built, the coordinates are transformed back to each of the 
distorted, cleaned FLT frames and PSF photometry is performed on each one of them. 
Positions and fluxes of each star are fitted simultaneously; also, stars with close neighbors 
are grouped together and fitted at the same time. Several passes are performed to eliminate 
the s mall contributions frora distant stars and to impr ove the init ial PSF model (derived 
from [Anderson fc Kine]|2004l with wings from Tiny Tim, lKristlll995l ) by calculating the PSF 
residuals. Aperture and charge-transfer efficiency corrections are applied to the flux from 
each individual FLT frame and the four values are finally combined to derive the observed 
VEGAMAG magnitude for each object. 

Several steps are needed to convert the observed magnitudes into masses. First, F550M 
values are transformed into the Johnson V band by adding 0.120 mag, as calculated for the 
appropriate input SED and extinction using the synthetic photometry package in CHORI- 
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Fig. 14.— ACS/HRC F550M DRZ frame of NGC 3603. The three WNha objects are marked 
and circles with radii of 3" and 6" centered on the cluster have been drawn. A square-root 
scale between and 800 counts has been used in order to show both the bright and the 
dim stars. The pixel size is 0'/025, the field size is 29'/3x25'.'4 (1172 px x 1016 px), and the 
vertical direction is 124° East of North. 
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ZOqlj. Second, Slogci — 5 is subtracted to obtain the absolute magnitudes V for each of 
the short, intermediate, and long distances. Finally, I use m{V) to obtain the corresponding 
apparent masses. The first three rows in Table 13.11 give the values I obtain for the three 
brightest objects in NGC 3603: Al, B, and C. 

In order to derive the AMF I first eliminate those objects with observed masses below 
8 M0. Using such a large cutoff allows us not to worry about incompleteness due to low 
S/N, since an 8 Mq isolated star has a S/N of 35 in a single FLT frame and of ~ 70 in 
the combined DRZ framj^. Also, since I want to compare the results with the experiments 
in the previous section, those stars with apparent masses above 120 (Al, B, and, for 
the two longest distances, C) are also excluded i.e. only the 8-120 Mq range is used for the 
power-law fit. No contamination correction is applied given the small size of the field, the 
high density of stars in the cluster, and that most such objects should be dimmer than the 
equivalent to the 8 Mq cutoff. Finally, the AMF is fitted using bins with an equal number 
of objects, as prescribed in Paper I. Different numbers of bins were tried but the slopes were 
left essentially unchanged, as expected from the results in that paper. The results for 10 bins 
are shown in Table 13.11 and plotted in Fig. [151 For comparison purposes, a fit with 10 bins 
of equal width was also attempted and, as expected, the obtained slopes were significantly 
different due to the bias inherent in such a method. Also, the measured uncertainties using 
bins of equal width were a factor of three larger (see Table 13.11) . Note that for the three 
distances 7equai number is Compatible (less than 1 sigma in two cases, less than 1.5 sigmas in 
the other) with the value we have used for our top-heavy experiments (7 = —2.0)0, so that 
will be the assumed slope for the rest of this subsection. 

Finally, to correct for the effects of unresolved binaries and/or chance superpositions, we 
need to calculate iVepx using Eqn. [TJ With the evolutionary synthesis module of CHORIZOS 
I obtain that the absolute V magnitude of a solar-metallicity 1-million-year old top-heavy 
"1-star cluster" is —0.75. The other quantity needed is the observed average magnitude in 
an effective pixel, which will depend on the aperture radius chosen. The DRZ field has a size 
of 29'.'3x25'.'4 but 69% of the flux is contained within a radius of 3" and 94% within a radius 
of 6". Therefore, the area enclosed within those two radii provide reasonable descriptions 



^''The value depends on the exact SED and Ay but the added uncertainty is of the order of only a few 
thousandths of a magnitude, which is good enough for our purposes. 

^■^Note that most incompleteness analyses deal with two different types of non-detections, those due to 
low S/N of the object and those due to the proximity of bright sources. Since the latter is already taken 
care of by our experiments, we only need to worry about the former. 

^^Of course, the observed 7 is that of the AMF, not the IMF, but, as we will see below, both should be 
very similar in this case. 
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Table 3. Distances and extinctions for NGC 3603. The uncertainties shown are standard 

deviations. 



Source d (kpc) Ay (mag) 



Moffat ri983) 


7.0 


± 


0.5 


4.46 


± 


0.28 


Melnick et ah ( 19891 


7.2 






4.47 


± 


0.40 


Eisenhauer et al. fl998) 








4.60 






de Free et al. f 19991 


6.1 


± 


0.6 








Niirnbereer et al. f2002l 


7.7 


± 


0.2 








Stolte et al. r2004l 


6.0 


± 


0.3 


4.50 


± 


0.60 


SuuEf & Bessell r2004l 


6.9 


± 


0.6 


4.44 






Melena et al. f2007l 


7.6 






4.66 







Table 4. Results for NGC 3603 for the three assumed distances. 



Quantity 


d = Q kpc 


d = 7 kpc 


d = 8 kpc 


Mai (Mo) 


159 


189 


217 


Mb (Mq) 


147 


176 


203 


Mc (Mo) 


105 


130 


152 


Tequal number 


-1.98 ± 0.16 


-1.98 ± 0.14 


-1.83 ± 0.12 


Tequal width 


-2.33 ± 0.49 


-2.68 ± 0.42 


-2.37 ± 0.39 


7lMF 


-1.90 ± 0.16 


-1.90 ± 0.14 


-1.74 ± 0.12 


iVepx (3") 


0.64 


0.88 


1.15 


iVepx (6") 


0.21 


0.29 


0.38 


^>12OM0,2s 


0.01 


0.01 


0.01 


^>12OM0,lb 


2.40 


2.83 


3.07 
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Fig. 15. — AMFs for NGC 3603 assuming a distance of 6.0 kpc (top), 7.0 kpc (middle), and 
8.0 kpc (bottom). In each case the data, the power-law fit in the 8-120 Mq range, and the 
top-heavy 2- binaries AMP are shown. Masses are in Mq. See the electronic version of the 
journal for a color version of this figure. 
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of the center of NGC 3603 and of the region that contains most of the hght from the core, 
respectively. The resulting A^epx are shown in Table 13. 1[ Depending on the distance and 
radius selected, I find that, on average, each point source in NGC 3603 contains at most 
1.15 individual stars due to chance superpositions with the assumption that all systems are 
single and not unresolved binaries. Therefore, chance superpositions are not an important 
issue for the NGC 3603 HRC data and, indeed, most effective pixels should be empty. Of 
course, if most point sources are binaries the results of the first experiment are still relevant. 

What do our experiments tell us about the corrections required to analyze the AMF 
in terms of the real IMF of NGC 3603? With the observed values of A^'epx, /a,ib (1 binary) 
should be the closest approximation to reality if indeed most stars are binaries (note that 
for binaries the values of iVepx have to be reduced by a factor of two and that those values 
represent the average number of binaries per effective pixel). As a comparison, I also discuss 
the case of fa,2s (2 single stars), which is less realistic but is the closest choice among the 
results of the second experiment. Regarding A7, both AMFs indicate that the measured 
slope for the range 8-120 will be slightly steeper than the real IMF (A7 = —0.086 fa,2s 
and by A7 = —0.087 for /a,ib)- Applying the correction — A7 we arrive at the final values 
for the IMF slope, 7imf; shown in Table 13. 1[ For the short and intermediate distances, I 
obtain a slope close to the —2.0 while for the long distan ce the slope is sligh tly steeper. The 



values are in good agreement with the recent result by IStolte et al.l (120061 ). who measured 
an IMF for NGC 3603 with 7 = -1.89 ± 0.14 between 7" and 20" (note, however, that their 
"Slower is significantly lower than 8 Mq). 

The result is different for the expected number of ultramassive objects, N^uoMq, shown 
in Table 13.11 For the single-star case, the expected number is almost zero. For the binary 
case, the predicted numbers are ~2 for the short distance and ~3 for the intermediate 
and long distances. In other words, if the real stellar upper mass limit is 120 Mq, we would 
expect two or three stars with apparent masses beyond that limit if most objects are unresolved 
binaries and that is precisely the number that is observed. 

W hat else do we know about Al, B, and C, the three brightest point sources in NGC 



3603? iDrissen et al.l (Il995l ) identified them as the source of the Wolf-Rayet part of the 



co mposite WR-l-0 integrated spectrum of the core, as had been already partially suggested 



by IWalbornI (119731 ). Since hydrogen is present and there is intrinsic absorption in the high- 
order Balmer lines, they are classified as WN6ha, implying that they are the more luminous 
and massive cousins of the 03 stars in their immediate vicinity. Therefore, from the point 
of view of their spectral classification they are good candidates to have masses above 120 
Mq if such stars really exist. However, it is also possible that they are unresolved multiple 
systems composed of individual stars with masses in the range 70-120 Mq which, from what 
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we cu rrently know, should have the same or similar spectral type (see e.g. iBonanos et al. 



2004 and note that their values are present-day masses, not initial ones). 

One way to decide which of the two opt ions is correc t is to search for radial velocity 



variations and possible eclipses, which is what iMoffat et al.l (120041 ) did. They found that Al 



is an eclipsing binary with a 3.7724 day period. Olivier Schnurr and Tony Moffat (private 
communication) have recently measured the radial velocity variations of Al and obtained 
masses of 114 ± 30 Mq and 84 ± 15 Mq, respectively. Interestingly, if we take those masses 
to be the initial ones, a 114 Mq + 84 unresolved object should have an apparent mass 
of 157 Mq] assuming some stellar mass loss during their lifetimes may add up to a few tens 
of solar masses to that value. The observed evolutionary masses lie between 159 Mq and 
217 Mq, which is in very good agreement with the Schnurr and Moffat result. 

Schnurr and Moffat (in preparation) have also found that C is a binary with an 8.92 
day period, so at least two of the WN objects in NGC 3603 are not single. B shows no 
signs of radial velocity variations so far but that does not mean that it is not a binary, 
since most massive binaries have separations (and some, also inclinations) that make their 
motions hard to de tect. On the other hand, both Al and C are overluminous in X-rays 



(iMoffat et al.ll2002l ). usually a sign of binarity since intense X-ray emission is produced in 
wind-wind interactions but B is normal. Then again, that does not preclude B from being 
a binary with a separation large enough not to emit extra X-rays but still small enough not 
to have its two components resolved by ACS/HRC. 

In summary, our analysis of NGC 3603 has found out that the existence of 2-3 objects 
with apparent masses above 120 Mq would be expected as the result of the observed number 
of stars in the cluster and its IMF below that limit without having to invoke the existence 
of stars with real masses above 120 Mq. The underlying hypothesis is that most or all 
massive stars in NGC 3603 are members of binary systems. Independent data have shown 
that indeed two of those stars are binaries and the currently available results for the third 
one leave room for the decision to go either way. The measurement of 7 could be revised 
slightly when the additional ACS data are processed. When we do that in a future paper, we 
will apply CHORIZOS to calculate individual extinction corrections for each star, derive a 
distance from the data, and do a comparison between the resulting color-magnitude diagram 
with the one in Fig. [H 
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3.2. R136 



From a historical point of view, R136 (HD 38268), the core of 30 Doradus, is an ex- 
cellent case stu dy on the proble ms caused by multiplicity, either real or induced by chance 



superpositions. IWalbornI (Il973l ). noting the WR features in its spectrum, compared it to 



NGC 3603 and suggested that it could be a Trap ezium-like system mad e out of Wolf-Rayet 
and other early-type O stars. On the other hand, ICassinelli et al.l (Il98ll ) dismissed that the 
central object, R136a, could be made out of "30 03 or WN3 within a space of 0.5 arcsec 
(0.1 pc)" and proposed instead that it was a single supermass ive star with 2500 M p^. A 



few years later that alternative was shown to be false fir st by iMoffat fc Seggewisd (Il983l ) 
by detecting spectral variations within R136 and later by IWeigelt fc Baied (119851 ) by using 
speckle interferometry to r esolve R136a into e ight individual sources, something which was 
later confirmed with HST (IWeigelt et al.lll99ll ). 



Several later HST st udies have analyzed the stellar composition of R136, of which the 
most complete is that of iMassey &: Hunterl (Il998l ). They obtained FOS spectroscopy of all 
the bright point sources in the core and combined it with PC (WFPC2) p hotometry to 



derive a color-magnitude diagram of the massive stars. In a follow-up paper, iMassey et al. 



( I2OO2I ) observed four spectroscopic eclipsing binaries and measured their masses. They also 
detected another five ec lipsing binaries, indicatin g that massive binaries are also common in 
R136. From their data, IMassey &: Hunteii (119981 ) measured a Salpeter- or Kroupa-like IMF 
with 7 between —2.3 and —2.4 for the range 15-120 Mq. That value, however, does not 
correspond to R136 alone but also includes the inner part of the 30 Doradus halo. In any 
case, for lack of better information, it will be the value assumed here. 

R136 is roughly seven times farther away than NGC 3603. Also, 30 Doradus has never 
been observed with ACS (either WFC or HRC) ar id the best data current ly available were 
obtained with the PC on WFPC2 (that was what IMassey fc Hunterl Il998l used), which has 
a factor of 1.6 lower pixel size than HRC. Therefore, if R136 and NGC 3603 had similar 
stellar densities at their centers, one would expect crowding to be (1.6 x 7)^ ~ 125 times 
worse in the former. The difference in resolution is readily apparent when one compares 
Figs. [H and [ 



Let us quantify how important crowding is within a l"-radius circle cen tered on R136. 
Doing aperture photometry on the F555W data, obtaining Ay from the data in IMassey &: Hunter 
( 119981 ) ■ and using CHORIZOS to estimate the conversion to the V passband, I derive an 
extinction-corrected mv of 9.29. Assuming an effective pixel radius of 2.5 real pixels I ob- 
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Fig. 16.— PC (WFPC2) F555W image of R136. The three WNha objects are marked and a 
circle with radius of 1" centered on the cluster has been drawn. A square-root scale between 
and 400 counts has been used in order to show both the bright and the dim stars. The 
pixel size is 0'.'0455, the field size is 4'.'23x3'.'64 (93 px x 80 px), and the vertical direction is 
37.8° West of North. The physical size (1.0 pc x 0.9 pc) is similar to that of Fig. [HI 
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tain mv,epx = 14.01. For assumed values of 1 million years (age)0, LMC (metallicity), and 
Kroupa (IMF), the value Myi = 0.79 is obtained. Finally, assuming a distance of 50 kpc, 
Eqn. [T] yields iVepx = 129. The value is relatively robust to small variations in the assumed 
parameters. Therefore, the n = 7 results of the second experiment (if binaries are not abun- 
dant in R136) and the n = 6 results of the third experiment (if binaries dominate) should be 
the reference for the analysis of the PC data of R136. This means that all the objects within 
the circle in Fig [16] have to be the superposition of a large number of stars (the majority of 
them being of low mass). 

Using /a,64b we obtain that 1.6% of the 77 effective pixels within 1" are expected to 
contain one object above 120 M©. The equivalent value for /a,i28s is 0.6%. This translates 
into 1.3 and 0.5 stars, respectively!^. How does this cor npare with the ob s ervati ons? Using 
the currently preferred low-temperature-scale values of iMassey &: Hunterl (119981 ) . the three 
most massive stars in the region are R136al, a2, and a3, with masses of 136 M©, 122 Mq, and 
120 Mq, respectively. Those observed masses agree very well with the values predicted by the 
experiments in this paper, especially if binaries dominate. It is also interesting to point out 
that three of the four mo st massive objects with spectral types in the larger region explored by 
Massey fc Hunteii (119981 ) are within the 1" radius considered here. One explanation would be 
that the most massive stars form preferentially in very dense environments. However, there 
are five other stars with similar spectral types (WNha and 03 If*/WN) outside the most 
crowded regions with masses only slightly lower (between 103 M0 and 120 M0). Hence, 
a reasonable alternative is that all of those stars have similar masses and that the extra 
measured mass in R136al, a2, and a3 is simply due to the increased likelihood of chance 
superpositions there. 



R136 was the cluster used by IWeidner fc Kroupal (|2004[) and iKoenI (120061 ) and one of 
the clusters in the sample used by lOey &: Clarkd (120051 ) to derive an mmax around 150 

(with varying degre e of uncertainty around the value). All of those works used the 
Massey &: Hunterl (119981 ) values for the masses which, even for the low-temperature scale, 
include at least a star with a mass of 136 M0 (R136al). Of course, the existence of a single 
star with such a mass invalidates the possibility that mmax is lower than that value. However, 



^^R136 could hardly be significantly younger than that because it has already cleared of dense gas a large 
region around it. Also, if its age were older than 2 million years, its most massive stars would no longer be 
of WNha type. 

^^Note that the expected number of ultramassive stars is lower for R136 than for NGC 3603, despite the 
facts that the regions considered here have the same extinction-corrected absolute magnitude and that R136 
is more crowded. The explanation resides in the different IMFs used for the calculation, Kroupa for R136 
and top-heavy for NGC 3603. 
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if, as we have seen, the certain existence of chance superpositions and the hkehhood that 
R136al is a binarjo (as most massive stars are) are considered, then the real mass of the 
heaviest component in R136al could be 120 Mq or even less. 

What about deriving the true IMF slope from the observed AMF at the center of R136? 
Figs. [8] and [13] indicate that A7 is in the range between 0.8 and 1.5 for 8-120 Mq in the 
region within 1". That correction is too large and too dependent on the assumptions to be 
applied with confidence, so with the current data it is not possible to know the real IMF 
slope at the very center of R136. However, note that if a factor of 2 in (linear) effective 
pixel size is achieved, a dramatic change takes place. The relevant functions then are /a,32s 
and /fe,i6s and A7 becomes ~ 0.10, which is a believable correction. Note also that in such 
a case, with ?a 300 effective pixels and with a S/N large enough to have 100-150 of them 
populated with detected stars, the number of observed objects will be enough to derive a 
bias-free meaningful IMF (a^ < 0.2) if the techniques in Paper I are used (see section 5 in 



Maiz Apellaniz fc Ubedal 120051 ). 



Our summary for R136 is similar to that for NGC 3603. Even though there are stars 
with apparent masses above 120 Mq, their detection can be explained by unresolved stars 
with real masses below that value. The main differences are that here the measured values 
are only slightly above the limit and that here the unresolved stars are a combination of 
chance superpositions and (likely) real multiple systems. 7 cannot be reliably measured 
with the current data for the innermost regions of R136 but a moderate improvement in the 
quality (within the capabilities of HRC/ACS) should be sufficient to do the job. 



4. Discussion 

I have presented in this paper a simple method to quantify the effect of unresolved 
multiple systems and chance alignments on the heavy part of the mass functions of young 
stellar clusters and applied it to NGC 3603 and R136. Unresolved multiple systems introduce 
an intrinsic spread in the observed color-magnitude diagrams and can produce a population 
of AUMSs in a rich cluster. On the other hand, their effect on the apparent massive star 
MF slope is relatively small. Chance alignments can also produce AUMSs but only for large 
Ns^p. As opposed to unresolved multiple systems, their effect on the AMF slope can be large. 

To my knowledge, no similar method had been previously developed for those same 
circumstances. For the effect of binaries in the IMF of low-mass stars in the field and in 



'Note that R136aH-a2 is not likely to be a bound system, as they are separated by at least 15 000 AU. 
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globular clusters, see iKroupa et al.l (Il99ll ) and ISollima et al.l (120071 ). respectively. For a 
similar (but more restricted) analy sis to the one in this paper applied to intermediate-age 
clusters, see lSagar fc Richtlerl (Il99ll ). Previous analyse s for young cluster s sometimes include 



a simplified correction with equal-mass binaries (e.g. iStolte et al.ll2006l ) but in most cases 
binaries are included in the incompleteness analysis. An incompleteness test adds a large 
number of fake stars to the observed field and runs the output through the star detection 
algorithm to see how many of the fake objects are recovered. The percentage of unrecovered 
objects is then used to calculate an incompleteness correction. Such an analysis in reality is 
measuring two different effects: the weakest objects (those close to the detection limit) may 
be missed, even in an uncrowded region, either because Poisson/background fluctuations 
or read noise place them below the detection limit or because the details of the detection 
algorithm make it stop before reaching it. Alternatively, some brighter objects are missed 
because they are located close to an even brighter star that hides them under its PSF. The 
flrst effect (low S/N objects) should be quantifiable by studying the detector properties and 
the expected S/N as a function of magnitude. The second effect is more complex because it 
depends on the PSF (which may be or may be not well characterized) and, more importantly, 
the stellar distribution. For that reason, incompleteness tests need to be run on a specific 
dataset. Furthermore, a classical incompleteness test does not address the issue that the 
added fake stars may significantly modify the measured magnitude of the detected, brighter 
star, and, thus, alter the results of the mass function. The technique presented here does 
not attempt to substitute the correction of the first effect (low S/N objects) provided by 
incompleteness tests but only the second. It does so by [a] correcting the measured slope 
after it is fitted to the observed mass function rather than modifying the observed mass 
function before fitting it and [b] including the effect of the different mass detected when 
two stars are superimposed and not only counting the fraction of undetected objects. It 
also has the advantage of being easily quantifiable and applicable to the planning of future 
observations. 

Having analyzed in detail NGC 3603 and R136, it is interesting to point out that no 
end-products of runaway collision s are seen there despite those two clusters being mentioned 
as likely candidates for them by iPortegies Zwart et al.l (120041 ) . There are several possible 
explanations: [a] Runaway collisions do take place but the remaining products evolve so 
fast that those in NGC 3603 and R136 have already become black holes (this would require 
lifetimes of ~ 1 million years), [b] Core collapse happens at a later stage, so the required 
time and density conditions are not satisfied in these clusters and ultramassive objects form 
only as a result of the merger of evolved objects, [c] Core collapse never happens and ultra- 
massive objects never form, [d] The output of a massive stellar merger that would lead to an 
ultramassive stars is so unstable that it loses its excess mass almost immediately. Regarding 
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the latter possibility, it is notorious that very massive stars are very difficult to model and 
that our current knowledge of them is limited. Nevertheless, we know that there could be 
different mechanisms at work that would impede the formation of an ultr amassive star via 



a me rger: the resulting object may be directly above Eddington' s limit (IG 



20071 ). so close to it that strange-mode instabilities overwhelm it (ITownsend 



always rotate so fast that it surpasses the QT limit (IMaeder fc Meynet 



ebbeek fc Pols 



200 



3 



20071 ). or may 



This article provides further evidence for the existence of a stellar upper mass limit at 
solar or near-solar metallicities. Furthermore, the analyzed data indicate that mmax may be 
as low as 120 Mq. The situation for Population III stars is expected to be different, likely 
because the red uced opacity should place Eddington's limit at higher luminosities and, hence. 



masses (see e.g. 



Zinnecker fc Yorkd 120071 and references therein) . 
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